# Load required packages and data
library(tidyverse)
library(pander)
library(fitdistrplus)
library(AER)
dat <- readRDS("t_test_data_sbs_pt01_10sub_10clust_10000sims.rds")

pval_b1_model <- unlist(dat[1,])
est_sb2 <- unlist(dat[2,]) # let's assume anything less than 1e-6 is the weird ones
tiny_est <- est_sb2 < 1e-6
resid_var_model <- unlist(dat[3,])
MSB <- unlist(dat[5,])
MSE <- unlist(dat[6,])
clustvar0 <- unlist(dat[7,])
clustvar1 <- unlist(dat[8,])
ftest_pval <- unlist(dat[9,])
pval_tt_eq_var <- unlist(dat[13,])
pval_tt_uneq_var <- unlist(dat[14,])
df_tt_uneq_var <- unlist(dat[15,])
beta0 <- unlist(dat[17,])
beta0_se <- unlist(dat[18,])
beta1 <- unlist(dat[19,])
beta1_se <- unlist(dat[20,])
plotdat <- cbind.data.frame(tiny_est, pval_b1_model, est_sb2, resid_var_model, #modelNLME,
                            MSB, MSE,
                            clustvar0, clustvar1, ftest_pval,# ttest_eq_var, ttest_uneq_var,
                            #clustmeans,
                            pval_tt_eq_var, pval_tt_uneq_var, df_tt_uneq_var,
                            beta0, beta0_se, beta1, beta1_se)#, test)
plotdat$icc <- MSB / (MSB + MSE)

Data-generating mechanism

Gathering data on how a ‘naive’ NLME model estimates those parameters, as well as what equal- and unequal-variance t-tests show.

When we get the tiny \(\sigma_b^2\) estimates, ratio of MSB to MSE from original dataset is small

Mean square within clusters: \[ MSE = \frac{1}{k(n-1)} \sum_i \sum_j (y_{ij} - \bar{y}_{i.} )^2 \] Mean square between clusters: \[ MSB = \frac{1}{k-1} \sum_i \sum_j (\bar{y}_{i.} - \bar{y}_{..} )^2 \]

I looked at a lot of relationships between variables, and the ratio \(MSB / MSE\) was the one that popped out as a clear predictor of a tiny estimate.

cutoff_MSBMSE <- min(MSB[tiny_est==F] / MSE[tiny_est==F])
cutoff_ICC <- min(plotdat$icc[tiny_est==F])

ggplot(data = plotdat) +
  geom_histogram(aes(x = MSB/MSE), color = "black", fill = "grey") +
  labs(title = "Ratio around .1 highly determinative of a tiny estimate") +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  geom_histogram(aes(x = MSB), color = "black", fill = "grey") +
  labs(title = "MSB alone is less predictive") +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  geom_point(aes(x = MSB, y = MSE, color = tiny_est), alpha = .2, stroke = 0) +
  labs(title = "Ratio around .1 highly determinative of a tiny estimate")

ggplot(data = plotdat) +
  geom_point(aes(x = MSB, y = MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_abline(aes(slope = 1/cutoff_MSBMSE, intercept = 0)) +
  labs(title = "More variability in that finding among the tiny estimates") +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  labs(title = "No change in the pattern as MSE varies") +
  geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_hline(aes(yintercept = cutoff_MSBMSE))

ggplot(data = plotdat) +
  labs(title = "No change in the pattern as MSB varies") +
  geom_point(aes(x = MSB, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) + geom_hline(aes(yintercept = cutoff_MSBMSE))

ggplot(data = plotdat) +
  geom_point(aes(x = MSB/MSE, y = clustvar0/clustvar1, color = tiny_est), alpha = .3, stroke = 0) + ylim(c(0,3)) +
  labs(title = "Pattern seems invariant to ratio of between-cluster variances")

ggplot(data = plotdat) +
  geom_point(aes(x = MSB/MSE, y = clustvar0 - clustvar1, color = tiny_est), alpha = .3, stroke = 0) +
  labs(title = "Pattern seems invariant to difference in between-cluster variances")

ggplot(data = plotdat) +
  geom_point(aes(x = MSE, y = icc, color = tiny_est), alpha = .2, stroke = 0) +
  labs(title = "Using the ICC instead of the ratio might also work, results are the same") +
  geom_hline(aes(yintercept = cutoff_ICC)) +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

ggplot(data = plotdat) +
  geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est), alpha = .2, stroke = 0) +
  labs(title = "Using the ratio") +
  geom_hline(aes(yintercept = cutoff_MSBMSE)) +
  facet_grid(rows = vars(tiny_est), labeller = label_both)

sum((MSB[tiny_est==T] / MSE[tiny_est==T]) > cutoff_MSBMSE)
## [1] 628
sum((plotdat$icc[tiny_est==T]) > cutoff_ICC)
## [1] 628
ggplot(data = plotdat) +
  geom_histogram(aes(x = pval_tt_uneq_var), color = "black", fill = "grey") +
  facet_grid(rows = vars(tiny_est), labeller = label_both) +
  xlab("p-value for welch t-test, allows unequal variances (results the same with eq vars)") +
  labs(title = "Small p-values on a t-test correlated with tiny estimates; makes sense, more power")

What to make of this? Just a signal-to-noise issue?

What about the cases that violate that boundary line?

First, find some wild outliers and examine them…

library(ggplot2)
outliers <- plotdat[ MSB/MSE > cutoff_MSBMSE+.05 & tiny_est == T, ]

outindx <- c(779, 4719, 6551, 7571)
plotdat$outlier <- FALSE
plotdat$outlier[outindx] <- TRUE

ggplot(data = plotdat) +
  geom_point(aes(x = MSE, y = MSB/MSE, color = tiny_est, stroke = outlier, alpha = outlier)) +
  labs(title = "Looking at these outliers") +
  geom_hline(aes(yintercept = cutoff_MSBMSE))
## Warning: Using alpha for a discrete variable is not advised.

outliers[,c(4:8,13,15)]
##      resid_var_model       MSB       MSE  clustvar0  clustvar1      beta0
## 779        0.9266005 0.1416865 0.9400950 0.04994531 0.10838582 -0.3350578
## 4719       0.9446727 0.1443756 0.9593768 0.03254286 0.12698338  0.1923321
## 6551       0.9413455 0.1562526 0.9545464 0.06520964 0.09665761 -0.3572094
## 7571       0.8345858 0.1554652 0.8560428 0.05251244 0.07149067  0.4466260
##           beta1
## 779   0.5034010
## 4719 -0.5113511
## 6551  0.5499081
## 7571 -0.6062691
ggplot(data = plotdat[tiny_est==T,]) +
  geom_point(aes(x = MSE, y = MSB/MSE, stroke = outlier, color = abs(beta1 - beta0) < .5), alpha = .5) +
  labs(title = "Outliers appear to show with large differences in model beta coefs") +
  geom_hline(aes(yintercept = cutoff_MSBMSE))

Something with the underlying data sets?

outlier_datasets <- dat[16,outindx]
ggplot(data = outlier_datasets[[1]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))

ggplot(data = outlier_datasets[[2]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))

ggplot(data = outlier_datasets[[3]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))

ggplot(data = outlier_datasets[[4]]) + geom_point(aes(x = clustid, y = linpred, color = arm_factor))

Trying A Tobit Fit

tdat <- cbind.data.frame(est_sb2)
tmod <- AER::tobit(est_sb2 ~ 1, left = 1e-6, data = tdat) # roughly log(-12)

# ggplot(data = as.data.frame(est_sb2)) +
#   geom_histogram(aes(x = est_sb2, y = ..density..), color = "black", fill = "grey", binwidth = .006) +
#   geom_vline(aes(xintercept = summary(tmod)$coefficients[1,1]), linetype = "dotdash", color = "red") +
#   geom_vline(aes(xintercept = sb2), color = "black") +
#   stat_function(fun = dnorm, args = list(mean = summary(tmod)$coefficients[1,1], sd = summary(tmod)$scale), color = "red") +
#   labs(title = "Tobit looks probable; note that mean of Tobit is lower than true value")

ggplot(data = as.data.frame(est_sb2)) +
  geom_histogram(aes(x = est_sb2, y = ..density..), color = "black", fill = "grey", binwidth = .006) +
  geom_vline(aes(xintercept = summary(tmod)$coefficients[1,1]), linetype = "dotdash", color = "red") +
  geom_vline(aes(xintercept = sb2), color = "black") +
  stat_function(fun = dnorm, args = list(mean = summary(tmod)$coefficients[1,1], sd = summary(tmod)$scale), color = "red") +
  ylim(c(0, 15))+
  labs(title = "Same graph, focusing on fit of non-zero parts")